#include "HHModel.h"
#include <cmath>


HHModel::HHModel():
        gK(36), gNa(120), gL(0.3),
        nuK(-12), nuNa(115), nuL(10.6),
        m(0.05), am(0), bm(0),
        n(0.3), an(0), bn(0),
        h(0.6), ah(0), bh(0),
        Iapp(0), Cm(1), nu(0), dt(0.1)
{
    setConstants();
}
void HHModel::setDt(double _dt)
{
    dt = _dt;
}
void HHModel::setIapp(double _Iapp)
{
    Iapp = _Iapp;
}
double HHModel::calculateNu()
{
    rungeNu();
    return nu;
}
double HHModel::getDt(void)
{
    return dt;
}
double HHModel::getM(void)
{
    return m;
}
double HHModel::getN(void)
{
    return n;
}
double HHModel::getH(void)
{
    return h;
}
double HHModel::getNu(void)
{
    return nu;
}

double HHModel::dm(double _m)
{
    return (am*(1-_m)-bm*_m)*dt;
}
void HHModel::rungeM()
{
    double k1,k2,k3,k4;
    k1=k2=k3=k4=0;

    k1 = dm(m);
    k2 = dm(m+k1/2);
    k3 = dm(m+k2/2);
    k4 = dm(m+k3);
    m = m + (k1+2*k2+2*k3+k4)/6;
}
double HHModel::dn(double _n)
{
    return (an*(1-_n)-bn*_n)*dt;
}
void HHModel::rungeN()
{
    double k1,k2,k3,k4;
    k1=k2=k3=k4=0;

    k1 = dn(n);
    k2 = dn(n+k1/2);
    k3 = dn(n+k2/2);
    k4 = dn(n+k3);
    n = n + (k1+2*k2+2*k3+k4)/6;
}
double HHModel::dh(double _h)
{
    return (ah*(1-_h)-bh*_h)*dt;
}
void HHModel::rungeH()
{
    double k1,k2,k3,k4;
    k1=k2=k3=k4=0;

    k1 = dh(h);
    k2 = dh(h+k1/2);
    k3 = dh(h+k2/2);
    k4 = dh(h+k3);
    h = h + (k1+2*k2+2*k3+k4)/6;
}
double HHModel::dnu(double _nu)
{
    return (-gK*pow(n,4)*(_nu-nuK) -gNa*pow(m,3)*h*(_nu-nuNa) -gL*(_nu - nuL) +Iapp)*dt/Cm;
}
void HHModel::rungeNu()
{
    setConstants();
    rungeM();
    rungeN();
    rungeH();

    double k1,k2,k3,k4;
    k1=k2=k3=k4=0;

    k1 = dnu(nu);
    k2 = dnu(nu+k1/2);
    k3 = dnu(nu+k2/2);
    k4 = dnu(nu+k3);
    nu = nu + (k1+2*k2+2*k3+k4)/6;
}
void HHModel::setConstants()
{
    am = 0.1*(25-nu)/(exp((25-nu)/10)-1);
    bm = 4*exp(-nu/18);

    an = 0.01*(10-nu)/(exp((10-nu)/10)-1);
    bn = 0.125*exp(-nu/80);

    ah = 0.07*exp(-nu/20);
    bh = 1/(exp((30-nu)/10)+1);
}
